Lobster Thermal Preference and GLORYS Bottom Temperature

# Libraries
library(tidyverse)
library(gmRi)
library(rnaturalearth)
library(scales)
library(sf)
library(gt)
library(patchwork)
library(ggdist)
library(ggimage)
library(tidyterra)
library(terra)
library(virtualspecies)
conflicted::conflict_prefer("select", "dplyr")
conflicted::conflict_prefer("filter", "dplyr")

# ggplot theme
theme_set(
  theme_gmri(
    axis.line.y = element_line(),
    axis.ticks.y = element_line(), 
    rect = element_rect(fill = "transparent", color = NA),
    panel.grid.major.y = element_blank(),
    strip.text.y = element_text(angle  = 0),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 8),
    legend.position = "bottom") +
    theme(
      legend.title = element_text(hjust = 0.5),
      legend.frame = element_rect(color = "black", fill = "transparent", linewidth = 0.8),
      legend.ticks = element_line(color = "black"),
      plot.background = element_rect(fill = "transparent", color = "Black"),
      panel.border = element_rect(fill = "transparent", color = "black", linewidth = 2), 
      legend.title.position = "top"
    ))



# Degree symbol
deg_c <- "\u00b0C"

# Gulf of Maine
gom_bounds <- read_sf(
  gmRi::get_timeseries_paths(
    region_group = "gmri_sst_focal_areas", 
   box_location = "cloudstorage")$apershing_gulf_of_maine$shape_path)


# Make a box to use when cropping based on an xlim and ylim pair
make_cropbox <- function(xlims, ylims){
  sfc <- st_sfc(st_polygon(list(
    rbind(c(xlims[[1]], ylims[[1]]),  
          c(xlims[[1]], ylims[[2]]), 
          c(xlims[[2]], ylims[[2]]), 
          c(xlims[[2]], ylims[[1]]), 
          c(xlims[[1]], ylims[[1]])))))
  sfc <- st_as_sf(sfc)
  return(sfc)
}

# Temp palette:
temp_pal <- rev(RColorBrewer::brewer.pal(n = 10, name = "RdBu"))


new_england <- ne_states("united states of america", returnclass = "sf")
canada <- ne_states("canada", returnclass = "sf")

# Use GMRI style
use_gmri_style_rmd()

Contents:

1. Progression of Temperature Change in NE Atlantic (annual, seasonal)
   
2. Representation of Warming Rates
   
3. Thermal habitat changes for lobster

Prepare Monthly Temperature Data

From a thermal preference angle, the suitability of certain habitats to lobster will change based on bottom temperature.

Now that we have access to freely available ocean reanalysis data it is possible to explore high resolution differences in bottom habitat suitability throughout the year.

Reference on Habitat Suitability

# Load GLORYS
gpath = cs_path("res", "GLORYs/NE_Shelf_TempSal_Monthly")

# Load the monthly averages and the climatology
# Done in: py/Monthly_Surface_and_Bottom.ipynb

# Temperature and anomalies
temps_monthly <- terra::rast(
  str_c(gpath, "Northwest_Atlantic_Surface_Bottom_93to2022_anoms.nc"))


# # Layer sources
# terra::sources(temps_monthly)
# terra::varnames(temps_monthly)


# Dates
r_dates <- terra::time(temps_monthly)

# # Why does it stack them this way?
# which(year(r_dates) == "2000" & month(r_dates) == 6)


# # climatology
# glorys_clim <- terra::rast(
#   str_c(gpath,"Northwest_Atlantic_Surface_Bottom_93to2022_clim.nc"))

# Define the bounding box as an sf object
# Cropping Extent  
crop_extent <- make_cropbox(
  xlims = c(-79, -56), 
  ylims = c(32, 51))

# Crop it
temps_cropped <- terra::crop(temps_monthly, crop_extent)
#varnames(temps_cropped)
bt_anom <- temps_cropped["bottom_temp_anom"]
sst_anom <- temps_cropped["surface_temp_anom"]
sst <- temps_cropped["surface_temp"][1:372]

# Split them into separate variables
# Find indices of layers that belong to "sst"
layer_names        <- names(temps_cropped)
sst_indices        <- which(str_detect(layer_names, "surface_temp"))
sst_anom_indices   <- which(str_detect(layer_names, "surface_temp_anom"))
bt_indices      <- which(str_detect(layer_names, "bottom_temp"))
bt_anom_indices <- which(str_detect(layer_names, "bottom_temp_anom"))

# and depth
depth_indices      <- which(str_detect(layer_names, "depth"))

# Subset the variables out separately by names
`%notin%` <- negate(`%in%`)
sst_indices <- sst_indices[which(sst_indices %notin% sst_anom_indices)]
bt_indices <- bt_indices[which(bt_indices %notin% bt_anom_indices)]

# Indexing is THE worst
sst      <- temps_cropped[[sst_indices]]
bt       <- temps_cropped[[bt_indices]]
sst_anom <- temps_cropped[[sst_anom_indices]]
bt_anom  <- temps_cropped[[bt_anom_indices]]
depths   <- temps_cropped["depth"][[1]]


# Assign dates as the names
names(sst) <- time(sst)
names(bt) <- time(bt)
names(sst_anom) <- time(sst_anom)
names(bt_anom) <- time(bt_anom)

# Take one depth layer to make mask for values over limit NA
bot_depth <- depths[[1]]
over_depths <- which(values(bot_depth) > 1450)
bot_depth[over_depths] <- NA

# # Plot the depth
# plot(bot_depth, main = "Max Depths")

# Decades
decades <- list(
  "avg_1990" = c(1990:1999),
  "avg_2000" = c(2000:2009),
  "avg_2010" = c(2010:2019)
)


# Get the mean temp for each decade
temp_decades <- purrr::map(decades, function(decade_yrs){
  
  layer_indices <- which(lubridate::year(time(sst)) %in% decade_yrs)
  
  sst_mean      <- terra::mean(sst[[layer_indices]], na.rm = T)
  sst_anom_mean <- terra::mean(sst_anom[[layer_indices]], na.rm = T)
  bt_mean       <- terra::mean(bt[[layer_indices]], na.rm = T)
  bt_anom_mean  <- terra::mean(bt_anom[[layer_indices]], na.rm = T)
  
  # Mask depths below our downloaded depth limit
  bt_mean[over_depths] <- NA
  bt_anom_mean[over_depths] <- NA

  decade_means <- list(
    "sst"      = sst_mean,
    "sst_anom" = sst_anom_mean,
    "bt"       = bt_mean,
    "bt_anom"  = bt_anom_mean)

  return(decade_means)
  
}) 

Decadal OISST

OISSTv2 Data is an observational dataset and not a reanalysis model and should be true-to-life.

# Load OISSt Monthly Means
sst_path <- cs_path("res", "OISST/oisst_mainstays/monthly_averages")

# Load the monthly averages

# Temperature and anomalies monthly
sst_monthly <- terra::rast(
  str_c(sst_path, "oisst_monthly.nc"))

layer_names <- names(sst_monthly)
both_indices <- which(str_detect(layer_names, "sst"))
only_anoms <- which(str_detect(layer_names, "sst_anom"))

# Subset the variables out separately by names
`%notin%` <- negate(`%in%`)
sst_indices <- both_indices[which(both_indices %notin% only_anoms)]

# I hate terra
sst_monthly <- sst_monthly[[sst_indices]]

# Dates
r_dates <- terra::time(sst_monthly)
names(sst_monthly) <- r_dates

# Decades
decades <- list(
  "avg_1990" = c(1990:1999),
  "avg_2000" = c(2000:2009),
  "avg_2010" = c(2010:2019)
)


# Get the mean temp for each decade
sst_decades <- purrr::map(decades, function(decade_yrs){
  
  layer_indices <- which(lubridate::year(time(sst_monthly)) %in% decade_yrs)
  
  sst_mean <- terra::mean(sst_monthly[[layer_indices]], na.rm = T)
  sst_mean <- terra::rotate(sst_mean)
  sst_mean <- terra::crop(sst_mean, crop_extent)
  decade_means <- sst_mean

  return(decade_means)
  
}) 

# Load ERSST
ersst_path <- cs_path("res", "ERSSTv5")
ersst <- terra::rast(str_c(ersst_path,"sst.mnmean_2022_02_02.nc"))
names(ersst) <- time(ersst)


# Do decadal averages for 1970's & 1980's
# Decades
decades <- list(
  "avg_1970" = c(1970:1979),
  "avg_1980" = c(1980:1989)
)


# Get the mean temp for each decade
ersst_decades <- purrr::map(decades, function(decade_yrs){
  
  layer_indices <- which(lubridate::year(time(ersst)) %in% decade_yrs)
  
  ersst_mean <- terra::mean(ersst[[layer_indices]], na.rm = T)
  ersst_mean <- terra::rotate(ersst_mean)
  ersst_mean <- terra::crop(ersst_mean, crop_extent)
  decade_means <- ersst_mean

  return(decade_means)
  
}) 

ERSSTv5 1970’s

ersst_17 <- ggplot() +
  geom_spatraster(data = ersst_decades$avg_1970) +
  geom_spatraster_contour_text(
    data = ersst_decades$avg_1970, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  map_theme(plot.margin = margin(1,1,1,1)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = T, crs = 4326) +
  labs(title = "Decadal Mean Surface Temperature\n1970's",
       fill = "Sea Surface Temperature")

ersst_17

ERSSTv5 1980’s

ersst_18 <- ggplot() +
  geom_spatraster(data = ersst_decades$avg_1980) +
  geom_spatraster_contour_text(
    data = ersst_decades$avg_1980, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  map_theme(plot.margin = margin(1,1,1,1)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = T, crs = 4326) +
  labs(title = "\n1980's",
       fill = "Sea Surface Temperature")

ersst_18

OISSTv2 1990’s

oisst_19 <- ggplot() +
  geom_spatraster(data = sst_decades$avg_1990) +
  geom_spatraster_contour_text(
    data = sst_decades$avg_1990, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  map_theme(plot.margin = margin(1,1,1,1)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = T, crs = 4326) +
  labs(title = "Decadal Mean Surface Temperature\n1990's",
       fill = "Sea Surface Temperature")

oisst_19

OISSTv2 2000’s

oisst_20 <- ggplot() +
  geom_spatraster(data = sst_decades$avg_2000) +
    geom_spatraster_contour_text(
    data = sst_decades$avg_2000, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  map_theme(plot.margin = margin(1,1,1,1)) +
 
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = T, crs = 4326) +
  labs(title = "\n2000's",
       fill = "Sea Surface Temperature")

oisst_20

OISSTv2 2010’s

oisst_21 <- ggplot() +
  geom_spatraster(data = sst_decades$avg_2010) +
  geom_spatraster_contour_text(
    data = sst_decades$avg_2010, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  map_theme(plot.margin = margin(1,1,1,1)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = T, crs = 4326) +
  labs(title = "\n2010's",
       fill = "Sea Surface Temperature")

oisst_21

Three Decades

three_decades <- (oisst_19 | oisst_20 | oisst_21) + plot_layout(guides = "collect") & guides(fill = guide_colorbar(barwidth = unit(6, "cm"))) 

three_decades

# # Save it
# # Prevents tiny fonts when saving
# showtext::showtext_opts(dpi=300) 
# 
# ggsave(
#   plot = three_decades, 
#   filename = here::here("local_data/images/OISSTv2_Decadal_Average_Maps_GMRI.png"),
#   height = unit(5, "in"),
#     width = unit(6, "in"),
#     dpi = "retina", 
#     bg = "white", 
#     scale = 1.75
# )

GLORYS Decadal Temperature Shifts

Changes in the decadal averages of surface and bottom temperatures for the region show a stark increase in temperature region wide in the 2010’s.

Anomaly values below are the departure from the 1993-2022 long-term average for each month, then averaged across the decades of interest.

SST

# Surface Temperature
sst_90 <- ggplot() +
  geom_spatraster(data = temp_decades$avg_1990$sst) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_1990$sst, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "1990's",
       fill = "Temperature")


# Surface Temperature
sst_20s <- ggplot() +
  geom_spatraster(data = temp_decades$avg_2000$sst) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2000$sst, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2000's",
       fill = "Temperature")
  


# Surface Temperature
sst_21s <- ggplot() +
  geom_spatraster(data = temp_decades$avg_2010$sst) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2010$sst, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(5,27)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2010's",
       fill = "Temperature")
  

(sst_90 | sst_20s | sst_21s) + plot_layout(guides = "collect") & guides(fill = guide_colorbar(barwidth = unit(6, "cm"))) 

SST Anomalies

# Surface Temperature Anomaly
sst_anom_90 <- ggplot() +
  geom_spatraster(
    data = temp_decades$avg_1990$sst_anom) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_1990$sst_anom,
    breaks = seq(-3,3,0.5),
    color = "gray40") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), oob = oob_squish) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "1990's",
       fill = "Temperature Anomaly")


# Surface Temperature Anomaly
sst_anom_20s <- ggplot() +
  geom_spatraster(
    data = temp_decades$avg_2000$sst_anom) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2000$sst_anom,
    breaks = seq(-3,3,0.5),
    color = "gray40") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), oob = oob_squish) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2000's",
       fill = "Temperature Anomaly")
  

# Surface Temperature Anomaly
sst_anom_21s <- ggplot() +
  geom_spatraster(
    data = temp_decades$avg_2010$sst_anom) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2010$sst_anom,
    breaks = seq(-3,3,0.5),
    color = "gray40") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), oob = oob_squish) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2010's",
       fill = "Temperature Anomaly")
  

(sst_anom_90 | sst_anom_20s | sst_anom_21s) + plot_layout(guides = "collect") & guides(fill = guide_colorbar(barwidth = unit(6, "cm"))) 

BT

# Bottom Temperature
bt_90 <- ggplot() +
  geom_spatraster(data = temp_decades$avg_1990$bt) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_1990$bt, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(0,16), na.value = "gray60", option = "magma") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "1990's",
       fill = "Temperature")


# Bottom Temperature
bt_20s <- ggplot() +
  geom_spatraster(data = temp_decades$avg_2000$bt) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2000$bt, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(0,16), na.value = "gray60", option = "magma") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2000's",
       fill = "Temperature")

# Bottom Temperature
bt_21s <- ggplot() +
  geom_spatraster(data = temp_decades$avg_2010$bt) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2010$bt, 
    breaks = seq(5, 30, 5),
    color = "white") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(limits = c(0,16), na.value = "gray60", option = "magma") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2010's",
       fill = "Temperature")



(bt_90 | bt_20s | bt_21s) + plot_layout(guides = "collect") & guides(fill = guide_colorbar(barwidth = unit(6, "cm"))) 

BT Anomalies

# Bottom Temperature Anomlay
bt_anom_90 <- ggplot() +
  geom_spatraster(data = temp_decades$avg_1990$bt_anom) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_1990$bt_anom,
    breaks = seq(-3,3,0.5),
    color = "gray60") +
  geom_spatraster_contour_text(
    data = temp_decades$avg_1990$bt_anom,
    breaks = seq(-3,3,0.5),
    color = "gray60") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), na.value = "gray60") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "1990's",
       fill = "Temperature Anomaly")


# Bottom Temperature Anomlay
bt_anom_20s <- ggplot() +
  geom_spatraster(data = temp_decades$avg_2000$bt_anom) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2000$bt_anom,
    breaks = seq(-3,3,0.5),
    color = "gray40") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), na.value = "gray60") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2000's",
       fill = "Temperature Anomaly")

# Bottom Temperature Anomlay
bt_anom_21s <- ggplot() +
  geom_spatraster(data = temp_decades$avg_2010$bt_anom) +
  geom_spatraster_contour_text(
    data = temp_decades$avg_2010$bt_anom,
    breaks = seq(-3,3,0.5),
    color = "gray40") +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_distiller(palette = "RdBu", limits = c(-1.5,1.5), na.value = "gray60") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  theme(panel.grid.major = element_line(color = "gray90")) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(title = "2010's",
       fill = "Temperature Anomaly")


(bt_anom_90 | bt_anom_20s | bt_anom_21s) + plot_layout(guides = "collect") & guides(fill = guide_colorbar(barwidth = unit(6, "cm"))) 

Lobster Thermal Preferences

Knowing that the temperatures globally are changing, it is logical to anticipate species to adjust their movements and behaviors to follow their thermal preferences.

From historical datasets and experiments we are able to determine the range of temperatures that species need and/or prefer in order to live.

Knowing these limits then lets us project onto a map where those conditions exist and where conditions are less-tolerable.

Lobster are a relatively well studied species, and it is believed that they prefer a range of temperatures between 10 and 20C. The effects of exposure to temperatures on either side of this range has an asymmetry such that lower temperatures are less of an issue than warmer ones. Below 10C it is understood that lobster activity and metabolism is lower, but they are otherwise unharmed. Above 20C temperatures have a stressful effect and can lead to disease and mortality.

# Get temperature range for study area as a vector
temp_range <- seq(
  min(values(temps_monthly["bottom_temp"]), na.rm = T),
  max(values(temps_monthly["surface_temp"]), na.rm = T),
  by = .1)


# Use the betaFun(), feed it temperature vector min/max and curve shape
lob_beta <- betaFun(
  x = temp_range,
        p1 = 0, 
        p2 = 20, 
        alpha = 4, 
        gamma = 0.7)


# Show the temp  preference
plot(
  lob_beta ~ temp_range, 
  type = "l", 
  main = "Hypothetical Lobster Temperature Preference Curve", 
  ylab = "Preference", 
  xlab = "Temperature Range")

Decadal Preference Maps

Lobster have a thermal preference (lab venture) of 11-22C. The following maps display what the decadal average temperatures look like with respect to the preference curve above.

# Run thermal preference for each decadal average


# Function to make preference rasters:
make_pref_ras <- function(in_ras, ras_name, p1, p2, alpha, gamma){
  
  # Make a new raster that we can swap values from  
  pref_ras <- in_ras[[ras_name]]

  # Assign values based on preference curve
  values(pref_ras) <- betaFun(
    values(pref_ras),
    p1 = p1,
    p2 = p2,
    alpha = alpha,
    gamma = gamma)
  
  return(pref_ras)

  
}




# Apply the preferences
lob_prefs <- map(temp_decades, function(decade_x){
  map(
    setNames(c("sst", "bt"), 
             c("sst", "bt")), 
    function(x){
      make_pref_ras(
        in_ras = decade_x, 
        ras_name =  x, 
        p1 = 0, 
        p2 = 20, 
        alpha = 4, 
        gamma = 0.7)
    })
  
})

Lobster 1990’s

sst_suit <- ggplot() +
  geom_spatraster(
    data = lob_prefs$avg_1990$sst) +
  geom_spatraster_contour(
    data = lob_prefs$avg_1990$sst, 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(
    title = "1990's",
    subtitle = "Surface Temperature Suitability",
    fill = "Temperature Suitability")


bt_suit <- ggplot() +
  geom_spatraster(
    data = lob_prefs$avg_1990$bt) +
  geom_spatraster_contour(
    data = lob_prefs$avg_1990$bt, 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  labs(
    title = "1990's",
    subtitle = "Bottom Temperature Suitability",
    fill = "Temperature Suitability")


(sst_suit | bt_suit)  + plot_layout(guides = "collect")

Lobster 2000’s

# Map of SST
sst_suit <- ggplot() +
  geom_spatraster(
    data = lob_prefs$avg_2000$sst) +
  geom_spatraster_contour(
    data = lob_prefs$avg_2000$sst, 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "2000's",
    subtitle = "Surface Temperature Suitability",
    fill = "Temperature Suitability")


bt_suit <- ggplot() +
  geom_spatraster(
    data = lob_prefs$avg_2000$bt) +
  geom_spatraster_contour(
    data = lob_prefs$avg_2000$bt, 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "2000's",
    subtitle = "Bottom Temperature Suitability",
    fill = "Temperature Suitability")


(sst_suit | bt_suit)  + plot_layout(guides = "collect")

Lobster 2010’s

# Map of SST
sst_suit <- ggplot() +
  geom_spatraster(
    data = lob_prefs$avg_2010$sst) +
  geom_spatraster_contour(
    data = lob_prefs$avg_2010$sst, 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "2010's",
    subtitle = "Surface Temperature Suitability",
    fill = "Temperature Suitability")


bt_suit <- ggplot() +
  geom_spatraster(
    data = lob_prefs$avg_2010$bt) +
  geom_spatraster_contour(
    data = lob_prefs$avg_2010$bt, 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75.8, -65), 
    ylim = c(36, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "2010's",
    subtitle = "Bottom Temperature Suitability",
    fill = "Temperature Suitability")


(sst_suit | bt_suit) + plot_layout(guides = "collect")

Monthly Thermal Preference Window

We can get a little more granular and apply that thermal preference window to monthly data. Below I’ve used 2023 as an example. This could also be done with daily data.

The following maps explore how different the thermal habitats for lobster can be from year-to-year by contrasting 2023 & 2024.

# Apply the activity logistic function
bot_suitability <- map(
  .x = setNames(names(bt), names(bt)),
  function(x){
      pref_ras <- make_pref_ras(
        in_ras = bt, 
        ras_name =  x, 
        p1 = 0, 
        p2 = 20, 
        alpha = 4, 
        gamma = 0.7)
      
      # Mask depths below our downloaded depth limit
      pref_ras[over_depths] <- NA
      return(pref_ras)
      
      
    })

May

map_23 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2023-05-31"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2023-05-31"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Temperature Suitability\nMay 2023",
    subtitle = "Lobster Suitability Range of 10-20C",
    fill = "Temperature Suitability")

map_24 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2024-05-31"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2024-05-31"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nMay 2024",
    subtitle = "",
    fill = "Temperature Suitability")

map_23 | map_24

June

map_23 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2023-06-30"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2023-06-30"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Temperature Suitability\nJune 2023",
    subtitle = "Lobster Suitability Range of 10-20C",
    fill = "Temperature Suitability")

map_24 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2024-06-30"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2024-06-30"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nJune 2024",
    subtitle = "",
    fill = "Temperature Suitability")

map_23 | map_24

July

map_23 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2023-07-31"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2023-07-31"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Temperature Suitability\nJuly 2023",
    subtitle = "Lobster Suitability Range of 10-20C",
    fill = "Temperature Suitability")

map_24 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2024-07-31"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2024-07-31"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nJuly 2024",
    subtitle = "",
    fill = "Temperature Suitability")

map_23 | map_24

August

map_23 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2023-08-31"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2023-08-31"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Temperature Suitability\nAugust 2023",
    subtitle = "Lobster Suitability Range of 10-20C",
    fill = "Temperature Suitability")

map_24 <- ggplot() +
  geom_spatraster(
    data = bot_suitability[["2024-08-31"]]) +
  geom_spatraster_contour(
    data = bot_suitability[["2024-08-31"]], 
    breaks = seq(0,1,.25), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nAugust 2024",
    subtitle = "",
    fill = "Temperature Suitability")

map_23 | map_24

Monthly Lobster Activity Level

Lobsters are thought to be relatively inactive below 10C. If we use this threshold we can map how activity appears spatially.

Use a logistic function centered on 10C, we can get a sense of how active lobsters are likely to be based on monthly bottom temperatures.

This could be adapted down the line to mark when* certain areas had their growing seasons begin and possibly some insight into their movements inshore/offshore.

For the below maps I am applying this logistic function to monthly bottom temperatures:

# Use the logisticFun(), feed it temperature vector min/max
lob_activity <- logisticFun(x = temp_range, alpha = -2, beta = 10)
plot(lob_activity ~ temp_range, type = "l", 
     main = "Lobster Minimum-Active Temperature", 
     ylab = "Activity Level", xlab = expression("Temperature"~degree~C))

# Function to make preference rasters:
logistic_pref_ras <- function(in_ras, ras_name, alpha = -2, beta = 10){
  
  # Make a new raster that we can swap values from  
  pref_ras <- in_ras[[ras_name]]

  # Assign values based on preference curve
  values(pref_ras) <- logisticFun(
    values(pref_ras),
    alpha = alpha, 
    beta = beta)
  
  return(pref_ras)

  
}





# Apply the activity logistic function
lob_bot_activity <- map(
  .x = setNames(names(bt), names(bt)),
  function(x){
      active_ras <- logistic_pref_ras(
        in_ras = bt, 
        ras_name =  x, 
        alpha = -2, 
        beta = 10)
      
      # Mask depths below our downloaded depth limit
      active_ras[over_depths] <- NA
      return(active_ras)
      
      
    })

May

map_23 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2023-05-31"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2023-05-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Activity\nMay 2023",
    subtitle = "Lobster Activity Beginning at Temps Above 10C",
    fill = "Lobster Activity Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2024-05-31"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2024-05-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nMay 2024",
    subtitle = "",
    fill = "Lobster Activity Level")


(map_23 | map_24) + plot_layout(guides = "collect")

June

map_23 <- ggplot() +
  geom_spatraster(
    data = lob_bot_activity[["2023-06-30"]]) +
  geom_spatraster_contour(
    data = lob_bot_activity[["2023-06-30"]], 
    breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Activity\nJune 2023",
    subtitle = "Lobster Activity Beginning at Temps Above 10C",
    fill = "Lobster Activity Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2024-06-30"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2024-06-30"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nJune 2024",
    subtitle = "",
    fill = "Lobster Activity Level")


(map_23 | map_24) + plot_layout(guides = "collect")

July

map_23 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2023-07-31"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2023-07-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Activity\nJuly 2023",
    subtitle = "Lobster Activity Beginning at Temps Above 10C",
    fill = "Lobster Activity Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2024-07-31"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2024-07-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nJuly 2024",
    subtitle = "",
    fill = "Lobster Activity Level")


(map_23 | map_24) + plot_layout(guides = "collect")

August

map_23 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2023-08-31"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2023-08-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Lobster Activity\nAugust 2023",
    subtitle = "Lobster Activity Beginning at Temps Above 10C",
    fill = "Lobster Activity Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_activity[["2024-08-31"]]) +
  geom_spatraster_contour(data = lob_bot_activity[["2024-08-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    labels = label_percent(),
    na.value = "gray60") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nAugust 2024",
    subtitle = "",
    fill = "Lobster Activity Level")


(map_23 | map_24) + plot_layout(guides = "collect")

Monthly Lobster Thermal Stress Maps

Lobsters are also thought to experience stress when under temperatures above 20C. At temperatures at or near these levels lobsters are likely to attempt to avoid these unfavorable conditions.

We can flip the logistic function and re-center it at/above 20C to then get a sense of where bottom temperatures may be too high.

For the below maps I am applying this logistic function to monthly bottom temperatures:

# Use the logisticFun(), feed it temperature vector min/max
lob_activity <- logisticFun(x = temp_range, alpha = -2, beta = 21)
plot(lob_activity ~ temp_range, type = "l", 
     main = "Lobster Thermal Stress Limit", 
     ylab = "Stress Level", xlab = expression("Temperature"~degree~C))

# Apply the activity logistic function
lob_bot_stress <- map(
  .x = setNames(names(bt), names(bt)),
  function(x){
      active_ras <- logistic_pref_ras(
        in_ras = bt, 
        ras_name =  x, 
        alpha = -2, 
        beta = 21)
      
      # Mask depths below our downloaded depth limit
      active_ras[over_depths] <- NA
      return(active_ras)
      
      
    })

May

map_23 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2023-05-31"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2023-05-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Thermal Stress Areas\n",
    subtitle = "Thermal Stress Beginning at Temps Above 20C",
    fill = "Stress Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2024-05-31"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2024-05-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nMay 2024",
    subtitle = "",
    fill = "Stress Level")

(map_23 | map_24) + plot_layout(guides = "collect")

June

map_23 <- ggplot() +
  geom_spatraster(
    data = lob_bot_stress[["2023-06-30"]]) +
  geom_spatraster_contour(
    data = lob_bot_stress[["2023-06-30"]], 
    breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
 guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
 labs(
    title = "Thermal Stress Areas\nJune 2023",
    subtitle = "Thermal Stress Beginning at Temps Above 20C",
    fill = "Stress Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2024-06-30"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2024-06-30"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nJune 2024",
    subtitle = "",
    fill = "Stress Level")

(map_23 | map_24) + plot_layout(guides = "collect")

July

map_23 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2023-07-31"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2023-07-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Thermal Stress Areas\nJuly 2023",
    subtitle = "Thermal Stress Beginning at Temps Above 20C",
    fill = "Stress Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2024-07-31"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2024-07-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nJuly 2024",
    subtitle = "",
    fill = "Stress Level")

(map_23 | map_24) + plot_layout(guides = "collect")

July

map_23 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2023-08-31"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2023-08-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "Thermal Stress Areas\nAugust 2023",
    subtitle = "Thermal Stress Beginning at Temps Above 20C",
    fill = "Stress Level")

map_24 <- ggplot() +
  geom_spatraster(data = lob_bot_stress[["2024-08-31"]]) +
  geom_spatraster_contour(data = lob_bot_stress[["2024-08-31"]], breaks = seq(0,1,.2), color = "white", linewidth = 0.15) +
  geom_sf(data = new_england) +
  geom_sf(data = canada) +
  scale_fill_viridis_c(
    breaks = seq(0,1, .2),
    na.value = "gray60",
    option = "plasma") +
  coord_sf(
    xlim = c(-75, -65), 
    ylim = c(39.5, 45.2), 
    expand = F, crs = st_crs(4269)) +
  guides(fill = guide_colorbar(barwidth = unit(4, "cm"))) +
  labs(
    title = "\nAugust 2024",
    subtitle = "",
    fill = "Stress Level")

(map_23 | map_24) + plot_layout(guides = "collect")